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§1 Billiards correlation functions. 



The first results on the ergodicity of the billiards go back about thirty years ago, [S]. But, as it is well 
known, ergodicity is a very weak, and in some sense, not too interesting, property. More direct physical 
meaning is attached to the correlation functions and to their decay speed. 

Let be the system phase space, St the evolution map, n(dx) the Liouville measure. 

In the case of the billiards is three dimensional: a point x £ is a point q £ M , where M is the billiards 
table, see Fig. 1, §3, i.e. a periodic box of side 1 with a few circles of radii Ri, R2, . . . taken out and regarded 
as obstacles, and an angle ip £ [0, 27r] so that (q, ip) represents the position and the velocity direction (with 
respect to the x-axis, say) of a point mass moving with unit velocity and direction ip until a collision takes 
place with the obstacles. Upon collision ip changes according to the elastic collision rules (equal incidence 
and reflection angles). The measure fi(dx) is simply fi(dx) = d 2 q d(p/(2ir\M\), where \M\ = area of M. The 
dynamical system (0, St, /J*) will be called the continous or 3d system (because it has a 3 dimensional phase 
space). The average with respect to [i over the phase space will be denoted by the symbol (.). 

An associated dynamical system is the collision system or 2d system: its (2 dimensional) phase space 
consists of points q £ dM and of the incidence angles d £ (f-, ^p) formed by the velocity at collision and 
the outer normal to the obstacle, counted counterclockwise. We denote (fi c , S c , /i c ) the "collision dynamical 
system" with fi c being the phase space, S c being the map mapping one collision £ = (s, d, i) to the previous 
collision £' = (s',d',i r ), where i is a label for the obstacle on which the collision takes place, s is the arc 
of the obstacle point q where the collision takes place (counted from an arbitrarily fixed origin), and d c = 
collision angle, and n c is the invariant measure: /i c (d£) = — cos d dd ds/normalization. The average with 
respect to /i c over the phase space fi c will be denoted by (.) c (and no confusion should arise). 

We shall consider three types of billiads (i.e. three configurations of obstacles): 00H, OH and D, defined 
precisely below. 

By correlations decay one does not mean the behaviour as t — > 00 of: 

</(<)/(0)> = J li(dx)f(Stx)f(x) (1.1) 

for the most general observable f, i.e. for the most general measurable function / on phase space. One is, 
in fact, interested in the (1.1) only for very special observables /. On the other hand it is clear that, the 
systems (0, St, /J*) and (fi c , S c , /i c ) being isomorphic to Bernoulli schemes, [GO], one can find nasty functions 
/ for which (1.1) approaches (/) as slowly as wished (by those who care about such wishes). 
A class (almost exhaustive) of interesting observables and related quantities is in the following list. 

1) the x component of the velocity v x = cos ip in the dynamical system (fi, S ( ,/i), leading to the velocity 
autocorrelation function: 

C(t) = (v x (t)v x (0)) (1.2) 
where v x (t) is the velocity at time t of an initial datum with velocity v x (0). 

2) the x component of the collision velocity, in the collision dynamical system (0 c ,5 c ,/i c ), v cx = cos ip where 
ip = d — ip q , ip q being the angle between the normal at the collision point q and the x-axis, (— 7r/2 < ip q < 
7r/2). This observable leads to the collision velocity autocorrelation function: 

Coin) = {v cx (n)v cx (Q)) c (1.3) 

where v cx (n) is the velocity at the n-th collision of an initial datum with velocity v cx (0). 

3) the transversal velocity autocorrelation given, with obvious notations, by: 

C T (t) = (v x (t)v y (0)) (1.4) 

Other related quantities are: 

4) the square displacements s(t) and s c (n), associated with the billiards and the collisions systems, respectively, 
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defined by: 

s(t) =(x 2 (t)) = 4 / dr f dT'(v x (T)v x (T')) 

Jo Jo ( L5 ) 
s(t) 

s c (n) ={x 2 (n)) c , D= lim 

t->co t 

where x(t), x(n) are two dimensional vectors which are measured by imagining to unfold on the plane the 
periodic box into all its images. And by thinking the motion to take place, without periodic boundary 
conditions but in the full plane, among the periodic lattice of scatterers generated by the obstacle images. 
In this case it is clear that D can be called the diffusion coefficient and the mean square displacement is, 
with the above notations, Dt (hence, for the purpose of comparison with [Bl] our diffusion coefficient D is 
related to the one, which we call D , in [Bl] eq. (1.9) by D = 2y/2D ). 

5) the number of collisions, v tm (x), undergone by the motion starting at x in the time t m . This is a random 
quantity if x is chosen randomly with distribution [i (or /i c ); then we can ask the value of the probability 
distribution: 

d a (z) = lim -L pro b ((i/ tm - (i/ t J)C e (z,z + 6z)) (1.6) 
t m ->oo dz I 

where {vt m ) is the average of v tm (x) over x and a > 0. 

The theoretical analysis of the above quantities started becoming at least conceivable after the work of [BS], 
about twelve years ago. 

The possibility of constructing a symbolic dynamics representation of the billiards motion, [BS], permitted 
to obtain the first upper bounds of the decay of C c (n). 

The bounds have been studied in the OH case (defined precisely below) and have the form: 

\C(t)\<ae- b(t ^°y, or |C c (n)| < e~ bn ' (1.7) 

for some e > (small), a, b > 0, where only the second bound is really proved in [BS]. 

The above (1.7) is just an upper bound and, very soon after [BS], the matter started being investigated 
numerically. The first results seemed to suggest that C(t) really did behave, on the sequence of times where 
the local maxima of |C(t)| are achieved, as a streched exponential with e ranging between 0.4 and 0.8 
depending on the scatterers geometry, [CCG], [BID]. A stretched exponential is a function of t or n of the 
form (1.7) and e is called the stretching. With the notable exception of [FM] who clearly stated that, in the 
OH case (see below), their results indicated an exponential decay of the sequence of local maxima, consistent 
with a representation of C(t) as a product of a pure exponential times a periodic function. 

Here by geometry of the scatterers one means some rough qualitative property: namely whether there are 
collisionless trajectories (billiards with horizon, denoted ooH), or not (billiards without horizon, denoted OH), 
when the periodic array consisting of the images of the obstacles does not allow to draw a path to infinity 
avoiding the obstacles, or diamond billiards, denoted D, in the other case where the obstacles keep the 
particle inside a bounded region. 

On the other hand the reason why one expects a stretched exponential decay seems to be related to the 
fact that the Markov partition realizing (see [BS]) the symbolic dynamical representation is not finite but 
countable. 

The denumerability of the Markov partition depends on the basic lack of smoothness of the billiards 
(0 c ,5 c ,/i c ) system dynamics: there are sharp discontinuities on the one dimensional lines corresponding 
to collisions which are preceded by a tangent collision (recall that we look at the motion backwards in time). 

If the Markov partition had been finite and the dynamical system smooth, the hyperbolicity of the bil- 
liards would have allowed us to conclude that the autocorrelation of any smooth observable / would have 
approached its n — > oo (infinite time) limit exponentially fast, in particular C c (n) would have been exponen- 
tially decaying (and hence presumably also C(t) because, under such (hypothetical) assumption, the time 
between collisions should be expected to be finite with finite moments). This is in fact a general mixing 
property for smooth observables in smooth hyperbolic systems, see [S2,Bow]. 

From the theory of the Markov partitions in [BS] one expects the partition to be "essentially finite", i.e. to 
consist of sets of apparently regular shape (irregularities breaking in two the elements of the partition, like 
tiny cuts, should be confined to very small scales), a finite number of which fills, for all practical purposes, 
the phase space. Furthermore, and most important, one expects to be able to construct "most" of the 
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elements of the partition by examining only the first few successive collisions of general initial data, at least 
in the cases of OH and D billiards, for which the time between collisions is finite with all its moments. 

If this is correct, as preliminary evidence deduced from our attempts to device a workable numerical al- 
gorithm to construct the Markov partition (in the OH case) seems to suggest, one would expect that the 
smooth observables like the velocity f = v x mix exponentially fast in the OH and D billiards: this means 
that in order to see non exponential behaviour (like stretched exponentials or even long time tails) one would 
need extremely accurate experiments over very many collisions. Otherwise the system should behave as a 
normal hyperbolic system without singularities. This is essentially our point of view (it could be called a 
"conjecture" but we refrain from formalizing it). 

In the case of ooH billiards, the same remarks should apply to the collision correlations, i.e. the correlations 
in the 2d collision dynamical system, even though the expected time between collisions has divergent moments 
of any order > 2. The correlations in the 3d dynamical system will have slow time decay (see [Bl]) as a 
consequence of the fact that there is also an infinite expectation value for the time of the first collision (as 
the initial velocities generating motions not experiencing collisions before a time t are those in contained an 
angle of order t -1 ). Note however that such a long time tail might become visible only for t so large to be 
not observable (with the present computer tools): and in fact this appears to be the case (see below). 

The above picture seems in sharp contrast with most of the existing numerical experiments, with the 
exception of [FM]. It does however agree with the theoretical work of [C] where a similar, but far simpler, 
non smooth hyperbolic system has been considered and exponential time decay of the correlations has been 
shown, in spite of the non existence of finite Markov partitions. The results of [C] led in fact Chernov to 
conjecture (independently) the exponential decay for the billiards collision systems as well. 

Having developed large sets of data, and seeing still quite far ahead the development of a numerical con- 
struction of a Markov partition, we thought that it would be useful to use our data to check the consequences 
of the above viewpoint on the correlation decay. Therefore we undertook a few extra experiments to see how 
strong was the evidence for stretched exponentials, and if one could confirm the important results of [FM], 
which had not received the deserved attention: unfortunately we became aware of this work too late for 
using exactly their triangular lattice configurations and to perform a detailed quantitative comparison. Our 
attempts at constructing Markov partitions had already led us to special billiard configurations which are 
not exactly the one of [FM]. Our work (being made about eight years after [FM]) contains better statistics 
(as computers are larger), but it is compatible with theirs. 

It should, however, be made clear, at the start, that in absence of rigorous asymptotics it is only possible to 
check if the experimental results are compatible with some a priori given function describing the decay. Thus 
it makes sense to test whether the experimental results can be fitted with a function of the form: e~ at f(wt) 
with / periodic. This implies that the successive maxima are placed on one or more parallel straight lines on 
a logarithmic plot. Of course one could test more complicated functional dependences (e.g. quasi periodic, 
aymptotically quasi periodic or stretched exponential times a periodic function, etc.), but the best one can 
expect from (our) analysis of the experimental results is a check of compatibility with an assumed law. If the 
assumed law allows for many parameters it is likely to give apparently better results; for instance a fit with: 
e -at f(wt) and / periodic, could give slightly better fit than the one with a = 1: this can be interpreted 
if, for instance, the best fit over (a, a) yields a close to 1, to mean that a is actually 1 and f(t) is close to 
periodic, with some small nonperiodic corrections. 

This difficulty with the analysis of experimental data is particularly evident in the case of the collision 
correlation function, C c (n), where the observations can only be made on integer values of the argument and 
even a periodic law would not necessarily yield data lying on parallel lines on a logarithmic plot. 

The following description of our experimental results seems to confirm the findings of [FM] and provides, in 
our opinion, further evidence for a purely exponential time decay of the correlation functions C(t), C c (n) as 
well as a diffusive behaviour of the motion (unfolded among the lattice of the images of the scatterers) when 
the scatteres do not confine the motion to a finite region (billiards ooH or OH, i.e. non diamond billiards). 
There are some puzzling features that we could not resolve, see comments to fig. 6. As a byproduct we 
studied also: 

1) The distribution of the number of collisions in a given time t m , for t m large, finding qualitative differences 
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between the ooH billiards and the OH or diamond billiards. 

2) The diffusion phenomena in the ooH and the OH cases. 

We shall (of course) argue that the "unexpected" exponential decay of the C(t) correlation function in the 
ooH case, and the associated diffusive behaviours is an artifact, because the long time tail is not yet visible 
on our computer experiment time scale. 

The diffusion seems "normal" (with mean square displacement of order t as t — > oo) in the OH case while 
is very likely (see [Bl]) anomalous (or "superdiffusive" with mean square displacement of order tint) in the 
ooH case: but such asymptotic regime does not seem to be visible at the time scales that we reach. 

Furthermore in the first case the distribution of the number of collisions occurring in a given time appears 
to be an asymptotically gaussian distribution. In the second case the distribution of the number of collisions 
appears to decay exponentially fast to the right and very sharply to the left of some value zo, although we 
again think that the latter is probably a short time effect. 

Our work seems to provide evidence for exponential decay of C(t) and C c (n): but the precision does not 
rule out exponential stretching with e close to one. A fortiori logarithmic stretching like: 

\C(t)\ < ae -K*/*o)(lnt/to)- s ( L8 ) 

is not ruled out either, and in fact it appears to us to be a natural candidate. But this can probably come 
out only from an analytic theory. 

The above discussion expresses numerical results obtained by experimenting over sequences of many col- 
lisions (more than the previous experiments) but still not in very large numbers (the hyperbolic nature 
of the motion prevents, even by performing double precision computations, to study more than about 20 
collisions in the D case, 15 in the ooH case and 10 collisions in the OH case, see fig. 5, for instance). What 
really happens beyond such time scales is not at present analyzable numerically, at least not in the sense of 
the present paper (and not much improvement can be expected in the future, even in the far future, as we 
are dealing with a chaotic system). Analytical work is here more promising, but of course difficult (see in 
particular [Bl]). 

One final comment: in order to avoid just repeating the work [FM] with better statistics we tried organize our 
data in such a way that they could be reliably reproduced (by those interested in them), including the error 
bars and the statistical analysis. This is explained in more detail in the next section, and it is the reason 
why we cannot consider very long time intervals. Remarks on numerical experiments. At the editors 
request we add a few general thoughts on numerical experiments on asymptotic properties of dynamical 
systems. There seems to be still concern on the actual measurability of quantities like Lyapunov exponents, 
decay of correlations, transport coefficients in numerical experiments. The problem seems to be that such 
quantities are defined by limits and "computers' cannot take limits". In fact the same argument could be 
construed for actual laboratory experiments: the latter seem to bother people less, as a true laboratory is 
far more reassuring than a computer, because it is older and we are used to it. It is important stressing that 
"computer experiments" are just experiments, worth of attention on their own; and the interpretation of 
their results may rest on theoretical frameworks which are more idealized. Measuring a Lyapunov exponent 
on a computer has the same meaning as measuring the period of the orbit of Jupiter, or the inclination 
of Mars. We measure them very accurately, but we are not sure that the motion is actually periodic, nor 
that the inclination is constant (they are not! but they "appear to be"). But in theoretical mechanics there 
are models of planetary motions in which Jupiter is moving periodically, or "can", and the Martian seasons 
are forever as dull as they are now. This gives us the "concept" of period and of inclination, and the idea 
of mesuring them, and of testing if they are constant. But of course our check of the constancy cannot 
be done as our lifes are too short (and the age of the universe is also too short): nevertheless we perform 
the measurements as accurately as possible and come out with precise figures. The figures tell us that, 
performing the measurements we did to measure the period or the inclination "as if they were well defined 
and constant", we get some definite results (by the way this may not be so easy). 

The interpretation of the results is customarily given in such a form: the period of Jupiter is T = ...±..% and 
the average inclination of Mars is i = ... ± ..%. Of course it is understood that they may change, or that they 
may simply not exist. But such doubt is not stated, usually; certainly not in the case or Mars' inclination. 
Hence this means that, on the basis of some theoretical models such motions are possible, and the data that 
we observe are consistent with the above values of T, i. Since the values are reproducible they qualify as 
"good data". What else could we possibly want? In the same sense in dynamical systems there are concepts 
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like mixing, ergodicity, Lyapunov exponents which are all asymptotic: they are well defined in abstract, 
idealized, models . We can however perform the operations that it would be necessary to measure them, 
get some results and then call such results with the corresponding names. We cannot check, in measuring 
Jupiter's period, that it returned after the jovian year exactly in the same position and not l.rara away 
(which could imply a double value of the period, and an infinite time to really check); likewise we do not 
check that the values of the Lyapunov exponents do not change by really going to a limit of infinite time or to 
twice the maximum time that we can reach (controlling the errors that we (think we) make). The numbers 
that we get are the "effective" Lyapunov exponents. In fact they are a far more interesting quantity than 
the abstract, non measurable, Lyapunov exponents: they are the ones "we would feel" if we had to interact 
with the system on the experimental time scale. The same can be said about the correlation functions. Of 
course the natural evolution of human affairs will perhaps make longer time scales experimentally reachable 
and new experiments to be performed. In fact it is now known that the average inclination of the axis of 
Mars is not fixed at all, it greatly moves, randomly (by about as much as 50°!), with a Lyapunov exponent 
of about the inverse of 1 million years, [LR]: which means that its motion might well be periodic (and with 
zero Lyapunov exponent) but for all practical purposes it is not such on rather short time scales (the result 
emerges from computer experiments, but once we finally succeed in being permitted to land on Mars, by the 
apparently unwilling inhabitants, we shall probably find geological evidence for that). 

In this paper we "measure correlations" in billiards. The deep theorems that are available do not give many 
hints, if any, on how to measure reliably (i.e. with an a priori controlled error) such quantities; in some cases 
there are not even theorems concerning their very existence, not to speak of approximability with finite 
algorithms (which is the goal, or should it be?, of reasonable theoreticians wishing to do scientific research 
and not philosophical speculations). Hence the real definition of the quantities is the one which emerges 
from the experiment that we describe. The minimum requirement is that it should be reproducible, and as 
accurate as possible: the language in which the results are formulated is necessarily borrowed from some 
general framework or "formalism", (ergodic theory and billiards theory in the present case). But they are 
not, and cannot, be theorems. They are just facts that hopefully, if the experiment is meaningful, may give 
new ideas to the theory and form a more satisfactory picture of the phenomena. But it is illusory that we can 
ever prove any theory by an experiment and viceversa: it is in fact incomprehensible to us why nevertheless 
a close connections between reality and theory can be at all empirically established. It might be because 
of what Galileo noted as the book of nature being written in mathematical characters: but books may not 
tell the truth, which in fact may just be undefined and undefinable. What really remains after performing 
an experiment is the expectation that others will find it interesting and useful (and this can only be if it is 
reproducible; and in spite of reproducibility it will often not be regarded as interesting or meaningful). 

We therefore use freely in what follows notions like correlations, distributions, decay, etc.: they are defined 
in the text, and ultimately in the present case by the computer program. And we tried to define clearly what 
we do, how we do estimate the errors; and we have avoided letting the computer just compute, trying to avoid 
the reader the frustation that we experienced many times in reading of numerical works (unquoted here) in 
which there was not enough information for reproducibility. We are, of course, aware that, for instance, the 
quantities like "diffusion coefficients", or "Lyapunov exponents", cannot in principle be measured from our 
data: but the values we provide for them are, for all practical purposes, the real ones in our experiment; 
and, if reproducible, in all experiments of the same kind. 

We know of no numerical experiments on billiards or other dynamical systems (even smooth) in which, for 
instance, the Lyapunov exponents are numerically estimated over a time longer than allowed to expand the 
round off error beyond the precision of a few percent (i.e. the usually reported errors on such exponents ): 
up to 20 collisions in our case. Therefore we have limited the duration of our experiments to a few collisions: 
more data could have been collected, but we could not have reliably intepreted them. For our billiards there 
are just not big enough computers for many more collisions (and there will never be, as their required size 
grows "exponentially" with the number of collisions). This is the reason it would be better to perform our 
measurements by other methods: but we know of no other methods for the billiards case. They might be 
developed in the future. We think that the pioneering age is over and the computer experiments should be 
subject to no less stringent standards than "ordinary" experiments (because we think they are no different). 
Therefore error bars should be mandatory, the error analysis should not be left out, round off should not 
be ignored, etc., and random number generators should be subject to tough tests. Many researchers do this 
and we would be happy to be considered at least close to them. 
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§2 The computer experiment. 



Our system is a square with periodic boundary conditions with sides of unit length, a = 1. We take the 
center of the torus as the origin of coordinates: (0, 0). There are a circle of radius Ri and center at (0, 0) and 
four more circles with radius R 2 and centers at (1/2, 1/2), (1/2, -1/2), (-1/2, -1/2), (-1/2, 1/2). Obviously 
only the part of the circles inside the torus is relevant (see Fig. 1). 




Fig. 1 

General billiard structure with scatterers of radius R\ and R2 in a box with side length a. 

A point particle is moving freely with unit velocity, \v\ = 1, in the space external to the circles and hitting 
them elastically (conserving the modulus of the total momentum and the energy). 
We studied three different cases (see Fig.l): 

Billiards with infinite horizon (00H): R\ = R2 = 7/20, less than the maximal value 4 _1 v / 2, but larger than 
the value 1/4 considered in [Bl], p. 366 (where the length unit is l/y/2, as the square lattice of the obstacles 
is with side 1 while in our case it is 1/V2; the velocity is however 1 in both cases so that, calling x(t) the 
position at time t, it is related to the position xsit) of [Bl] by x(t) = -^XB(tV2)). 

Billiards without horizon (OH): R x = 1/5, R 2 = 2/5. 
Diamond (£>): R 1 = R 2 = (5/32) 1 / 2 > 4" 1 v / 2. 

The value for the radius chosen for the Diamond case is the same of [CCG] and it will make it easier the 
comparison with their results. 

We use the following algorithm. The initial particle position is chosen at random with uniform distribution 
from the external space left by the circles. We give to the particle an initial velocity vector with an angle with 
respect to the horizontal axis chosen from a uniform distribution between and 2w. We compute whether 
the future particle trajectory will hit a circle or a box side and, in every case, the particle is moved up to 
one of those possibilities. If the particle hits a box side, periodic boundary conditions are applied and if it 
hits a circle we change the velocity direction according to the equations for an elastic collision with a surface 
with a given curvature. We apply the algorithm until the time of the particle evolution is equal to 4 units 
of time (an average of 20 collisions, as we shall see) and then we start again the algorithm. 

We do not apply molecular dynamics: the succesive collisions are determined by solving, on the computer, 
the appropriate equations for the geometric intersections. 

Note that for the OH and D billiards, the above experiment permits us to study the collision correlations 
without performing new experiments. In fact, by choosing initial 3d data with distribution [i automatically 
produces a distribution ji c of the first collisions. And we performed a few experiments measuring the collision 
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(2d) correlation function even though they are very difficult to interpret (as commented in §1). 

For the ooH case, in contrast, by choosing initial 3d data it is not possible to produce the distribution 
H c because some collisions require to occur a time larger than our experiment. Therefore, in the case of 
ooH, we have also done a set of computer simulations starting with an initial datum which is a collision 
(i.e. the position is randomly choosen on an obstacle boundary with uniform distribution and its direction 
is a collision angle d chosen randomly between 7r/2 and 37r/2 with distribution — cos d dd). Hence in the 
latter case we study both, the 2d system and the 3d system. This is necessary because, as said above, for 
the ooH case the average distance between collisions in the 3d system is oo, while in the 2d case it is finite 
(although it has a divergent second moment). Therefore, a priori, the latter property could affect the time 
behavior of the macroscopic averages. On the other hand, it is very useful in itself to study the 2d case in 
order to compare with the most recent theoretical works on this billiard type (p.e. see [Bl]). 

In all the graphs referring to the 3d systems the time is measured in units of: 

to = I<;f d8dc ( 2 .!) 
J dsdc 

where dc = — cos d dd and t(s, $) is the time between the collision (s, $) and the previous one; so that to 
is the theoretical average time between collisions (note that, by the ergodicity theorem [S], this is exactly 
equal to the average time between collisions in the 2d collision system as dsdc is the invariant measure in 
the collision system defined at the begining of §1). In particular, we get: to = 0.151..., 0.316... and 0.1645... 
for D, OH and ooH cases. 

The plotted functions were the average over 2 x 10 7 trajectories for the Diamond case, 1 x 10 7 trajectories 
for ooH case and up to 2.8 x 10 8 trajectories for the OH case. This implies a statistical error proportional to 

N -l/2 _ 2 . 1Q -4 for D 3 . 1Q -4 for ^jj and g . 1Q -5 for qjj 

In order to provide reproducible data and errors we describe in the Appendix our precise definition of the 
type of fit we do on the data, the evaluation of the errors on the fitting parameters and the criteria used 
to measure quantitatively how good different fits are. We think that our data are careful enough to be 
reproducible, errors included. 

We were worried about the rounding error propagation due to the system chaoticity. This error is propor- 
tional to s(0)e At where s(0) is the initial rounding error (10 -16 , in our case) and A is the system maximum 
Lyapunov exponent. We know from our own computer simulation estimates that A ~ l.OStQ 1 for D , 
A ~ 1.75*0 1 for 0H and A - 1-23 1^ 1 for ooH. The original results existing in the literature are, in the D 
case, due to [Be] and our computations agree within the errors with them. 

Then, we may conclude that the statistical error dominates over the rounding error when t, measured in 
absolute units of time, is less than 4.2, 5.2 and 3.9 for D, OH and ooH cases respectively. Therefore, in our 
computer simulation (going up to a time (measured in absolute units) < 4), the only relevant source of errors 
should be the statistical ones. Computer simulations with longer time interval will measure hydrodynamical 
long time properties, i.e. time correlations between "different" initial and final trajectories. 

The error bars reported in the graphs are the mean square dispersions of each reported value added to an 
estimate of the other error sources (round off propagated by chaoticity). They are always plotted although 
most of the times they cannot be seen because of the graph width. 

Since we are interested in asymptotic quantities we do not want to include in our fits the short time data 
as they certainly show transient phenomena. Therefore we decided, quite arbitrarily, to discard the data 
produced in the first two units of time measured in units of the appropiate Lyapunov exponents. Such time 
scale, Lyapunov time scale, is different from the time scale used in the graphs (namely the (2.1)): therefore 
we marked in each graph the Lyapunov time scale by an arrow pointing at its value in the x-axis. 

In the 3d system (i.e. continuous time dynamical system) we measured the magnitudes at times t, = 
i * t ra /4000 (i = 1, 2, .., 4000), where t m is the measured time interval in absolute units. The computation of 
such 3d magnitudes is the most expensive part in CPU computer time because for each trajectory and variable 
we have to perform, at least, 4000 operations more than the usual ones for the 2d system (approximately 
1500 operations in a trajectory of 15 collisions). We represent in the correponding graphs 400 points that we 
got by averaging locally 10 data points. The averaging is done only to reduce the number of data plotted and 
it does not differ appreciably from the full plot (i.e. the first it is within the error bars of the second) which 
represents the data really used in our fits. The symbols regularly used in the plots are small black circles (or 
ellipses) and big empty ones for the 3d and 2d functions respectively. The error bars appear (when visible) 
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as vertical bars. 



§5 Billiards velocity autocorrelations: C(t) = (v x (0)v x (t)} . 



a: The ooH (3d) case: 



In Fig. 2 we show the |C(t)| behavior. We see the characteristic oscillations with monotonically decreasing 
amplitude. We manage to observe up to 6 or 7 oscillations before the statistical errors obscure the data. We 
realized how the amplitudes decrease with an apparently regular law and the oscillation period seems to be 
constant r = 2.55(±0.08) t - 

1. I ' 1 ' 1 ' 1 ' 1 



\C(t)\ 



0.8 r 



0.6 - 



0.4 - 



0.2 - 



0. 




12. 



16. 

(ooH) 



Fig. 2 



Absolute value of \C(t)\, the velocity-velocity correlation function, versus t/to for the ooH case. The arrow marks the Lyapunov 
time scale (here and in the following pictures). 



Hence we try the fit of the data with a law e~ at l to f(wt) with / 27r-periodic. If this is a good guess we 
should be able to get a good fit of the data In |C(<) | , evaluated on the relative maxima points t,, with a 
function of the form: —at /to + b, with b = hif(wti). If the function / is not too complicated, the time t, 
should appear at multiples of 2irw~ 1 . More generally one could expect the maxima to occur on parallel lines 
of the form —at/to + with 6, = ln/(ij) and t{ separated by 2irw~ 1 . 



This is confirmed when we plot the In |C(i)| versus t (see Fig. 3). We discard the points with t < Zto k, 2A _1 
and we may fit a straight line crossing the last 5 maxima: — 0.34(±0.05) — 0.345(±0.006)t/to This implies 
a pure exponential law behavior for the amplitude with a decay rate 2.90(±0.05)to> apparently unrelated to 
the Lyapunov exponent. 
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0. 

-2. 
-4. 

ln|C(<)| 

-6. 
-8. 
-10. 

0. 4. 8. 12. 16. 

t/t (ooH) 

Fig. 3 

Logarithm, of the velocity-velocity correlation function absolute value, ln|C(t)|, versus t/to for the ooH case. The solid line is 
the best linear fit of the last 5 local maxima of \C(t)\. The fit starts at twice the Lyapunov time scale (see comments at the 
end of $2). 

We also tried a three parameters best fit with the stretched exponential fit for the maxima getting an 
exponent 0.93 with error bars which does not include the unit value. Of course the fit is somewhat better 
and this is not surprising (clearly had we allowed for two more free parameters we would have obtained a 
perfect fit, see remarks in §1). 

b): The OH case: 

We performed an analysis similar to the above. In Fig. 4 we see that |C(i)| has more structure than in the 
ooH case due to the more complex structure of the unit cell of the lattice of obstacles, (consisting of two 
sublattices with circles with different radius). 

1. 
0.8 
0.6 

\C(t)\ 

0.4 
0.2 
0. 

0. 4. 8. 12. 

t/to (OH) 

Fig. 4 

Absolute value of \C(t)\, the velocity-velocity correlation function, versus t/to f or the OH case. 

The maxima of |C(i)| have an apparently non trivial structure and this suggest a e~ at l to f(wt) form of the 
curve with / 27r-periodic but more complicated in structure. This can be tested by a (rather severe) two 
parameter fit and the maxima should lie on parallel straight lines of the form at, /to + hi in a logarithmic 
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plot. 

When we plot In |C(i)| versus t (see Fig. 5), we see that there is in fact such a regular behaviour: there are 
two well defined similar cells with three maxima in each one. We get the fits for the two pairs of contiguous 
maxima: -0.9(±0.3) - 0.54(±0.05)i/*o and -1.2(±0.1) - 0.54(±0.04)i/*o- In both cases the decay rate for 
the correlation function is 1.8(±0.1)io> apparently uncorrelated to the Lyapunov exponent. Note that the 
error bars, at the maxima points, become too large after t/to ps 12. 




Fig. 5 




This experiment is very time consuming: we performed it with great care over a much longer time span 
than the others. The reason is that an experiment over a time scale comparable (computerwise) to the others 
would have yielded Fig. 5, which would have been quite inconclusive because of the error bars becoming 
too large after 9to (with the corresponding data not used in the fit). On the other hand the large number of 
relative maxima promised a very accurate test of our exponential decay assumption if the calculation could 
be pushed to the time 12io> where in Fig. 5 the linear law looks possible only because of the large error 
bars. The longer time experiment results allows us to reduce error bars considerably up to t ~ 12io- We 
reported the results also for t larger that 12io> when the error bars become visible just to give an idea of 
how difficult it could be to improve the analysis (which in the present form took three months of CPU on 
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an IBM RISC-6000, just for Fig. 6). 

Although the more refined experiment confirms and improves all the data of Fig. 5, the last clear maximum 
in Fig. 6 seems to be much higher than it should be in order to fit on the lower line. This might be due 
to corrections to the periodicity hypothesis or to actual stretching of the decay; or our definitions of errors 
are too generous (because, see Appendix, they (must) contain some elements of arbitrariness, as in (A1.6), 
which in such extreme data may require further analysis). This is a point that imposes further care, as we 
could not resolve it. 
c): The Diamond case: 

In Fig. 7 we show |C(i)| in this case. Its structure is similar to the ooH case: oscillations with decreasing 
amplitude with period r = 1.72(±0.03)^o- In this case we have very clean (errorwise) data up to times of 
order 13 t/to as it can be seen more clearly from the following logarithmic plot. 

1. 



0.8 - 



\C(t)\ 



0.6 " 



0.4 - 



0.2 - 



0. 




Fig. 7 

Absolute value of \C(t)\, the velocity-velocity correlation function, versus t/to for the D case. 



\n\C(t)\ 




Fig. 8 

Logarithm, of the velocity-velocity correlation function absolute value, ln|C(t)|, versus t/to f or the D case. The solid line is 
the best linear fit of the last 8 local maxima of \C(t)\. 

The plot of In |C(i)| versus t shows again the pure exponential behaviour of the last 8 maxima amplitudes 
that we fit by the equation (Fig.8): 0.28(±0.06) - 0.49(±0.01)i/*o- The decay rate for the correlation 
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function is now 2.04(±0.04)^o, again apparently unrelated to the Lyapunov exponent. As in the ooH case, 
if we try to fit a stretched exponential, we get an exponent 0.995 with an error which does not include the 
value 1 (see comments in §1 and to Fig. 4 above). 
Comments: 

There are not too many results in the literature for this autocorrelation function. In the ooH case [BID] 
get, from a computer simulation with a single trajectory with 2 • 10 7 collisions, that: |C(i)| ~ e~ at for 
t < to, and |C(i)| ~ e~ bt for to < t < 43^o- and, furthermore, referring to the decay, found in [FM]: 
|C(i)| ~ ct~ x when t m 110. to they remark that it is independent of R when R < 2 -3 / 2 . They also give 
arguments for a pure exponential decay when the system dimension goes to infinity. There is no reference 
to the rich structure of C(t), clear in [FM], in the analysis of the decay. 

The paper [FM] also contains, as mentioned in the introduction a large number of important results, that 
our experiments confirm; we mention: 1) They found for an "almost" triangular diamond a behavior similar 
to ours in Fig. 8 and they said explicitly that: "the velocity autocorrelation function appears to be an 
exponential damped periodic function". 2) They showed for the OH case a rich structure for maxima and 
minina less ordered but similar to ours in Fig. 5 and 6. They say explicitly that "the curve appears to be 
bounded by an exponential envelope". 3) They did not show for the ooH case the regular decay as we show 
in Fig. 3 for short times. But they show a rough consistency of a power law decay as l/t for long times. 

§^ Collision velocity autocorrelation: C c (n) = (v cx (0)v cx (n)) c , 

The results of this section are a byproduct of the numerical results obtained by performing the experiments 
described in the previous section with the exception of the ooH case. Although we do not think that they are 
particularly significant for the reasons explained in §1, we report them as their study illustrates the problems 
on the data analysis on the collision systems stressed in §1. 

a): The ooH (3d) case: 

We computed the velocity autocorrelation function corresponding to the hit number n. For large n the noise 
due to statistical errors is important. Since the computer simulation was performed for a fixed time t m = 4. 
(in absolute units), each trajectory may have different number N(n) of collisions for each n in time t m ; the 
number of data N(n) that we averaged to get C c (n), depends on n and so does its statistical error (because 
for large n there are large fluctuations in the number of the sampled trajectories that actually experienced 
n collisions, see also §7). We find that for n > 16 collisions the errors (which, as always, have two natures: 
statistical and dynamical) become too large to permit setting up a data analysis. 

0. 

-2. 

-4. 

ln|C c (n)| 

-6. 

-8. 

-10. 

0. 6. 12. 18. 

71 (ooH) 

Fig. 9 

ln|C c (n)| versus n in the ooH (3d) case. The straight lines are the best linear fits for the data points crossed by the lines. 
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The Fig. 9 shows a In |C c (n)| versus n plot. We see that there is a regular behaviour for the even and odd n 
values separately. The two sets of values may be separately fitted by pure exponential laws as it is shown in 
the figure. The fits are: -1.28(±0.03) - 0.228(±0.004) n and -2.22(±0.06) - 0.171(±0.007) n for the upper 
and lower lines respectively. The decay rate for each set is 4.28(±0.08) and 5.8(±0.2) respectively. Also it 
seems that one of the sets merges in the other when n > 15. 

The latter singular behavior may be perhaps understood by the discrete nature of C c (n). For instance, 
if the continuous time function had the analytical form: C(t) = e~ at f (wt / to) , being / a periodic function 
with period w~Ho, (see §3a) then it would seem natural to expect that the collision one has the form: 
C c (n) = e~ an f(w'n). Assuming that w' is near a "resonance", say wo = h/k, with h,k integers, and 
expanding the function / around wo we get: In |C c (n)| = —a'n + In f(won) + (w' — wo)nf (won)/ f(won). 
Therefore, we can get k straight lines as the function /'// nas period k, at least if (w' — wo)n is small. 

In [BID] C c (n) is computed by using a unique trajectory with 2-10 7 collisions (we computed 2-10 7 trajectories 
with 30 collisions each (at least), i.e. 6 • 10 8 colllisions). They fitted C c (n) = (—l) n e~ a " h for 1 < n < 9 with 
b = 0.86 ± 0.06. They do not mention the two branches that we observe in C c (n). We do not understand 
well enough their results to be able to reproduce or interpret them. 

The OH case: 

In Fig. 10 we show the ln|C c (n)| versus n behaviour. In this case there do not seem to exist different 
branches for n even and odd as they appeared in ooH (it could be due to the fact that the period of 
oscillations may be far from a resonant regime and therefore, in the observed scale, the phenomenon is not 
relevant enough to be appreciated). All points fit the line: 1.1(±0.3) — 1.2(±0.1)n (see Fig. 10). 




Fig. 10 

In |C c (n)| versus n in the OH case. The line is the best linear fit for the data points it touches. 

Note how in this case the interval that can be analyzed only goes up to n ps 10 compared with n ps 20 
for the ooH case. That is so because the mean free time in this case is about twice the one for the ooH 
case (because of our choice of the geometrical parameters) and therefore the number of collisions with good 
statistics is reduced by a factor of two. Even though in Fig. 10 it is not shown, the data follow the rule: 
sign(C c (n)) = (-1)" ■ 

The D case: 

The results in this case are less clear. In Fig. 11 we show ln|C c (n)| versus n. The data, compared with 
the previous cases, have fluctuating behavior. We get some linear fits by using different n intervals with 
relative errors in the coefficients larger than 15%. For instance, the one showed in Fig. 11 is: — 1.(±2.) — 
0.5(±0.3)n. We also get a three parameters stretched exponential fit (see comments in §1 and following 
Fig. 2) starting the minimization process randomly (with the parameters values ranging between —5. and 5.) 
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we found: 4.0(±0.9) — 4.1(±0.4)n°' 41 ( ±ao3 ). In order to compare with the non stretched exponential law 
we fix the value of the stretched exponent to 0.41 and we fit the remaining two parameter function. We 
get: 4. (±2.) — 4.(±l.)n ' 41 . The ratio between the goodness (see Appendix) of the exponential and stretched 
exponential behaviors is: 1.1. Plotting in a 3d graph the error function G which measures the accuracy of the 
fits (see Appendix), as a function of the three parameters of the stretched exponential fit, say (ao> <*i> <*2)> 
being «2 the stretched exponent, one finds that there is a path connecting the values (4.0,-4.1,0.41) and 
(— 1., —0.5, 1.) which is almost flat, slowly going up from G = 0.4281 to G = 0.4688. We conclude that the 
avalaible data are by far not good enough to distinguish stretched exponential laws with stretching parameter 
«2 anywhere between 0.4 and 1. 

0. 



-4. 

ln|C c (n)| 

-8. 



-12. 

0. 6. 12. 18. 

(D) 

Fig. 11 

In |C c (n)| versus n in the D case. The line is the best linear fit for the data points lying in the same n-interval of the line. 

In [CCG] this case is studied by averaging 2.5 x 10 5 different trajectories and computing the quantity 
r(n) = (£a(0)£a(^)) where ^(n) is the characteristic function of set A at time n and A is the set: 

A = {(x,y,v x ,v y \xe [0,1/3], y £ [0,1/2], cos(l.l) <v x < cos(O.l) 
and sin(O.l) < v y < sin(l.l)} 

where the arguments of the trigonometric functions are measured in radians. The data are fit to the stretched 
exponential r(n) = e~ 1An . Two main points from their computer simulation are not clear enough to us: 
there is an oscillatory character in their data which they eliminate by using a "suitable smoothing procedure", 
i.e. they seem to eliminate the "structure" of the data (a structure that might require further analysis as 
Fig. 3 and Fig. 8 illustrate). On the other hand they study trajectories with 70 collisions, and we cannot (as 
commented above) go beyond about 20 collisions still hoping to control the rounding error propagation and 
the statistical error which is too large to have good data. It would be necessary to know better the numerical 
method and data analysis actually used in the experiment [CCG] to clarify the above matters. 



§5 The transverse autocorrelation C T (t) = (^(O)^^)}. 



We also computed the transverse velocity self correlation function for the ooH (3d) case. In Fig. 12 and 
Fig. 13 we plotted its time behaviour: absolute value and the logarithm the absolute value versus time. The 
amplitude of oscillations is two order of magnitude smaller compared with the one in C(t) and their period 
is equal. Again, the maxima follow an exponential decay behavior. In Fig. 13 we see again how the maxima 
decay in a regular fashion and the linear fit of the last four maxima is: — 2.27(±0.07) — 0.37(±0.02)i/io- The 
decay rate is equal to the one of C(t) within error bars, i.e. 2.70(±0.15)io- 
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We have also computed the C T (t) in the OH and D cases, but we do not report the data (as no new 
information comes from them). 



The diffusion coefficients and s(t) = (x 2 (t)), s c (n) = (x 2 (n)) c . 



a) The ooH (2d and 3d) case. 

In Fig. 14 we show the mean square particle displacement, s(t), versus time for the 2d and 3d systems. In 
both cases and for short times there is some structure due to the correlation with the initial condition (in 
particular we see the expected parabolic behavior coming from the free motion when the time is near zero). 
Notice also how the 2d system has a larger mean square displacement because its mean path to the first 
collision is larger. The latter explanation looks like inconsistent with the existence of unbounded trajectories 
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for the 3d system. In fact, there exist trajectories with free paths larger than the time interval used 
(i.e. t m ) but their number is so small that they, effectively, don't contribute to the averages. The asymptotic 
lines for the 2d and 3d cases are respectively: 0.069(±0.04) + 0.0077(±0.0003)i/*o and 0.0369(±0.0006) + 
0.00758(±0.00004)V*o- 

0.20 I ' 1 ' 1 ' 1 ' 1 



0.15 " 



s(t) 0.10 



0.05 - 




Fig. 14 

Mean square displacement, s(t), versus time t, for the ooH (3d and 2d) case. The equations are the best linear fits for the 2d 
and 3d data which is larger than 14. to 
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Fig. 15 

Mean square displacement, s c {n), versus the collision number, n, for the ooH 3d and 2d cases. The equations are the best 
linear fits for the 2d and 3d data which is between 8 and 16. 

Recently [Bl] argued that for the 2d system, the mean square displacement has the superdiffusive behavior: 
s(t) = a + -Dstln(t) + Dt when t is large enough. We tried to fit this behavior to our data but we found 
that D B ~ 0.0001*0 1 and D ~ 0.0075* o ~\ where the error is in the fourth significative digit (i.e. 100% !). 
That is, the time interval we are using is too small to detect clearly the asymptotic superdiffusive behavior. 
In fact, we are far from the asymptotic regime in order to compare with the theoretical results by [Bl]. The 
latter 
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remark is corroborated when we get Db following its analytical expression given in [Bl]: 



Db = -^-D° B (V2R) (6.1) 

and: 

D B (x) = ^-(l-2xf (6.2) 

having converted the results to our units (with lattice spacing 1) and R is the obstacle radius (in the same 
units). In our case it gives the value Db = 0.00014^ 1 which is of similar order of magnitude of our computer 
simulated value. 

The Fig. 15 also shows s c (n) versus n. The asymptotic lines for the 2d and 3d cases are respectively: 
0.097(±.007) + 0.0063(±.0004)n and 0.0720(±.0007) + 0.00591(±.00005)n. If we fit the funtion s c (n) = 
a' + D' B n\n(n) + D'n + b'n~ x to the data we find the same result for the a! and D' parameters, being 
D' B zero within our error bars. Finally, the relation D/D' gives: 1.22(±0.13)i o and 1.28(±0.011)i o - The 
discrepancy, far beyond the mean square deviation, with respect to the expected result s(t)/s c (n) = to when 
t, n go to oo ([Bl]) possibly indicates again that we are far from the asymptotic region. 

The OH case: 

In Fig. 16 and Fig. 17 we show s(t) and s c (n) respectively. We get in both cases asymptotic linear 
behaviours: a + Dt/t with a = -0.026(±0.003), D = 0.1129(±0.0004)*o 1 and a' = 0.027(±0.001), 
D' = 0.1072(±0.0002) and then D/D' = 1.053(±0.006)i o - We see that the error bars does not include 
the expected exact result D/D' = to: we think that is because we are not yet in the asymptotic regime. 

Note we do not fit the last three data points in Fig. 17. This is so because the time interval used (t m = 4) 
is relatively short and some particle trajectories do move long enough to experience more than 8 collisions 
(i.e. they are trajectories with long free paths) and therefore they are relevant for the average but they are 
not contributing to the final value of s c (n) (in other words, at absolute time t m = 4 there are trajectories 
which have not yet experienced n collisions; their number is negligible if n < 9 but starts becoming important 
if n > 9, see Fig. 21 below). 
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Fig. 16 

Mean square displacement, s(t), versus time t, for the OH case. The equation is the best linear fit for the data which is larger 
than 6 to 
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s c (n) 




Fig. 17 

Mean square displacement, s c {n), versus the collision number, n, for the OH case. The equation is the best linear fit for the 
data with n is between 4 and 9. 



c): The D case: 



In Fig. 18 and Fig. 19 we show s(t) and s c (n) respectively. We get the expected non-diffusive behavior: 
s(t) and s c (n) tend to the limiting constant values 0.017 and 0.026 respectively. We can compare, to check 
the ergodicity, such numerical values with the ergodic averages: 



f dxd-d x 2 f dsdc x 2 

s(+ ^ = L JdZdT = °- 01727 -' s ^ = L Jd^r 



0.0264. 
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Fig. 18 

Mean square displacement, s(t), versus time t, for the D case. The number is the asymptotic exact value (see text). 
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acement, s c {n), versus the collision number, n, for the D case. The number is the asymptotic exact value 



This is a nice check of the ergodic theorem ([S]) and of course it could be repeated for many other averages 
and used as a test of how much one could expect to be in an asymptotic regime. Such a test would probably 
be necessary if one wanted to perform a more complete analysis of the ooH and OH cases. In the D case, as 
we see, there are no reasons to expect that the asymptotic regime has not yet been reached in the observed 
time intervals. 



§7 The collision number distribution: d(n;t m ). 



We computed the probability distribution, d(n;t m ), that the particle undergoes n collisions in the time 
interval [0,i ra ] for t m = 1, 1.5, 2, 2.5, 3, 3.5 and 4. In Fig. 20, 21 and 22 we show d(n;t m ) versus n for the 
ooH (3d), OH and D cases respectively. 
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Fig. 20 

The distribution probability of the number of trajectories having n collisions in a time interval [0,t m 
(3d) case for (from left to right) t m = 1., 1.5 ? 2., 2.5 ? 3 V 3.5 and 4.; the dotted lines are visual guides. 



d(n;t m ), in the ooH 
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The aim of this analysis was to understand the reliability of the calculations on the ooH case. We also 
performed the same calculation for the OH and D cases (see below); in the latter cases the computation was 
done in order to compare with the ooH case. As a byproduct we could test in the three cases the natural 
hypothesis that a suitable rescaling leads to a gaussian distribution of the number of collisions in a given 
time. 
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Fig. 22 

The distribution probability of the number of trajectories having n collisions in a time interval [Q,t m ], 
t m in the OH (Fig, 21) and D (Fig. 22) cases for (from left to right) t 




d(n',t m ), for different 



We see that, except in the ooH (3d) case, the distributions look like gaussians and, in each case, their 
structure evolves with t m in a regular form. In particular in Fig. 23 we show the averaged number of collisions 
per trajectory, (n) ( , for all three cases and different t m . We see how the linear fit is almost perfect in all 
cases: 0.9998(±0.0001) + 6.622(±0.004)* ra , 0.99836(±0. 00008) + 6.081(±0.004)* ra and 1.0005(±0.0004) + 
3.228(±0.002)t ra for the D, OH and ooH cases respectively. The inverse of the linear fit slope gives in all 
cases the mean free time (within errors) to for each case respectively. 
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Fig. 23 

The average number of collisions by trajectory, (n) f , in a time interval [0,t m ] for the ooH (3d), OH and D cases. The dashed 
lines are the best linear fits of the data points. 

In Fig. 24 we show the logarithm of the distribution standard deviation (as computed from the data in 
the graphs in Fig. 19), \na tm , for all three cases and different t m . The fit ailn(t ra ) + 02 gives: (01,02) = 
(0.4(±0.1),-0.04(±0.06)) for™ the D case, (0.51 (±0.04), 0.35 (±0.02)) for the 00H case and (0.51 (±0.04), 
0.35 (±0.02)) for the OH case. 
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Fig. 24 

The logarithm of the standard deviation of the probability distribution d(n;t m ) versus ln(t m ). The dashed lines are the best 
linear fit of the data points. 

From the above scaling information, it is natural to expect that the distribution of the random variable: 
z = (n — (n) ( )/(T( m defined by d a (z) = a tm d(n(z);t m ) will be independent of t m . 

This is shown in Fig. 25, Fig. 26 and Fig. 27 where we plotted all data for different t m in the 00H (3d), OH 
and D cases respectively. 

In the 00H (3d) case, the scaled distribution, d a (z), is not symmetric around the maximum. The distribution 
right wing (i.e. from the maximum to n — > 00) seems to have a pure exponential behaviour. The left wing 
has a very fast decay. But it seems that the distribution depends on t m indicating that we are far from the 
asymptotic regime (which might still be a gaussian). For the OH and D cases the scaled data distribution 
can be fitted by a gaussian distribution: as expected from the rigorous theory, see [BS], [BSC]. 
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Fig. 27 

-\- <7t m z; t m ) with different t m for the D case, (z) is the average values of 



§ Appendix. Fits and Errors. The random number generator. 



From our computer experiments we get data sets corresponding to averages over many dynamical trajecto- 
ries, say y(x) = {j/(£i)}i=i,jv, where x = {x,}, = i t N is in our case an independent variable, say for instance 
a set of TV" collision numbers or time instants. To the latter experimental data set of points, we want to 
fit a given guessed function, say f(x; a), where a = {a„} n= i tP is a set of arbitrary parameters. Here by fit 
we mean to find a set of parameters a* which optimizes some reasonable functional relation between the 
experimental data and the fitting function. 

In our case we use the least squares funtional, i.e. : 

N 

V(y(x), a) = J2 [Vi ~ /(^! «)] 2 ( A1A ) 

i = l 

The set of parameters a*(y) is here obtained by asking that they should be the minima of the V function: 
ds*V(y, a*) = . We also define the goodness of our fit, G, by the average j/-distance of our data to the 
function f(x;a*): G(y(x)) = V(y(x), a*)/N. This parameter is only meaningful when it is compared with 
the one from another fit. Given many fits, the one with smallest G value will be called best fit (among the 
considered fits). 

The data have, in general, non-negligible errors, say e = {Ei}i = i t N, due to the finite number of samples 
used in the averaging and to the rounding error propagation because of the system chaoticity (see comments 
in §2). Such errors induce errors on the parameter values. Therefore, a measure of the error amplitude in 
a*(y) is given by: 

A e ~a*(y) = [a*(y + e)-a*(y-e)}/2 (A1.2) 
In the particular case in which the magnitude of the data error is much smaller than the measured value, 
\^i/y( x i)\ << 1) we may expand the latter equation around e = 0: 

N 

A e ^„(# = $>(">(£)<■,■ (41.3) 

! = 1 

The coefficients cf 1 ^ are found by expanding V(y(x) + e, a) around a*(y) and e = and they are given by: 

c\ n) = -J2 D^ n d 2 y(xi)a ,J(y, a*) (A1A) 

m = 1 
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where D mn = <9^ a , V(y, a*). 

In particular for the linear fit, f(x, a) = a± + «2 x, the coefficients c^'^ 2 ^ are given by: 



where Ax = (x — x) 2 and x n = TV" -1 X^=i K i 
The errors are random variables and we have to average them over their distribution. Since our data we 
have come from dynamical trajectories, the errors may be correlated and we cannot assume that they are 
independent gaussian distributed random variables. Therefore we empirically estimate an upper bound for 
their correlation values: 

< ^(£l)(£ 2 j)e~ lxi ~ XjlX (41.6) 

where A is the corresponding Lyapunov exponent. The parameter errors in our analysis are defined by the 
equations: 



N N 



i=i i=i 

where 6f = (e?) and their use and meaning is described entirely by the above comments. 
Finally a comment on our random number generator: we have used the so called R250 in which a sequence of 
pseudo-random numbers, {X(n)}, is generated by the linear recursion relation: X(n) = X(n— 103). xor.X(n— 
250) where .xor. is the logical operation exclusive or along the bits of both numbers, i.e. l.xor.l = 0, 
O.xor.O = and l.xor.O = 1. The first 250 random numbers are generated with a random number modulo 
generator: X(n) = 16807 *X(n- l)mod(2 31 - 1). The recurrence period for the R250 random number is 
expected to be (2 103 - 1) * (2 250 - 1). We do not have reasons to think that the random numbers could be 
"bad"; our results should be reproducible, although their interpretation might be different from ours. 
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